Generate data with imperfections

Fit normal HMM with 2 states

Try normal model with 3 states

model <- stan_model("hmm.stan")

data_stan <- list(
  N=nrow(df),
  dist=df$obs,
  K=3,
  cost=10,
  p_anomalous=0.01,
  conc_anomalous=100
)

fit <- optimise_repeat(5, model, data_stan)

df$state_est <- fit$par$state
df$state_est[1] <- 1

g <- df %>%
  mutate(state_est=as.factor(state_est)) %>% 
  ggplot(aes(x=time, y=obs)) +
  geom_line() +
  geom_point(aes(colour=state_est)) +
  xlab("Time") +
  ylab(latex2exp::TeX("$\\Y_t")) +
  scale_color_brewer("State", palette = "Dark2")

g


# ggsave("../figures/simulated_corrupt_normal_state_3.pdf", g, width = 6,
#        height = 3)

Fit a model with rubbish collection state.

data_stan <- list(
  N=nrow(df),
  dist=df$obs,
  K=2,
  cost=10,
  p_anomalous=0.01,
  conc_anomalous=100
)
model <- stan_model("hmm_bernoulli.stan")
fit <- optimise_repeat(5, model, data_stan)

df$state_est <- fit$par$state
df$state_est[1] <- 1

g <- df %>%
  mutate(state_est=as.factor(state_est)) %>% 
  ggplot(aes(x=time, y=obs)) +
  geom_line() +
  geom_point(aes(colour=state_est)) +
  xlab("Time") +
  ylab(latex2exp::TeX("$\\Y_t")) +
  scale_color_brewer("State", palette = "Dark2")

g


# ggsave("../figures/simulated_corrupt_normal_state_3.pdf", g, width = 6,
#        height = 3)
fit$par$phi
[1] 0.01184509
fit$par$mu
[1] 0.9560906 8.0321222
fit$par$sigma
[1] 1.198069 1.321352

2D HMM

Fit 2D robust HMM

data_stan <- list(
  N=nrow(df),
  dist=x %>% as.matrix(),
  K=2,
  cost=7,
  p_anomalous=0.1,
  conc_anomalous=100
)
model <- stan_model("hmm_bernoulli_2d.stan")
fit <- optimise_repeat(10, model, data_stan)

g <- x %>% 
  mutate(state=fit$par$state) %>% 
  ggplot(aes(x=V1, y=V2)) +
  geom_point(aes(colour=as.factor(state))) +
  xlab(latex2exp::TeX("$\\Y_{1t}")) +
  ylab(latex2exp::TeX("$\\Y_{2t}")) +
  scale_color_brewer("State", palette = "Dark2")
g


# ggsave("../figures/simulated_2d_corrupt_rubbish.pdf", g, width = 6,
#        height = 3)
fit_r <- fit
fit$par$mu
            [,1]        [,2]
[1,]  7.88259312  7.96074249
[2,] -0.04429909 -0.08280655
fit$par$sigma
[1] 1.015293 1.038413

Try standard HMM with 2 states

fit_2 <- fit
fit$par$mu
           [,1]       [,2]
[1,]  7.5813273  7.6541182
[2,] -0.2256035 -0.2759647
fit$par$sigma
[1] 1.541746 1.581182

Try standard HMM with 3 states

fit_3 <- fit
fit$par$mu
           [,1]      [,2]
[1,] -0.3135811 -0.372200
[2,]  3.5225315  3.660930
[3,]  7.8946867  7.964446
fit$par$sigma
[1] 1.295070 1.335942

Plot state cleanliness

mu_r <- fit_r$par$mu
mu_r <- mu_r[c(2, 1), ]
mu_2 <- fit_2$par$mu
mu_2 <- mu_2[c(2, 1), ]
mu_3 <- fit_3$par$mu[c(1, 3), ]
mu_true <- matrix(c(0, 0, 8, 8), ncol=2, byrow = TRUE)

tmp <- tribble(
  ~y1, ~y2, ~state, ~label,
  mu_true[1, 1], mu_true[1, 2], 1, "truth",
  mu_true[2, 1], mu_true[2, 2], 2, "truth",
  mu_r[1, 1], mu_r[1, 2], 1, "rubbish",
  mu_r[2, 1], mu_r[2, 2], 2, "rubbish",
  mu_2[1, 1], mu_2[1, 2], 1, "2-state normal",
  mu_2[2, 1], mu_2[2, 2], 2, "2-state normal",
  mu_3[1, 1], mu_3[1, 2], 1, "3-state normal",
  mu_3[2, 1], mu_3[2, 2], 2, "3-state normal"
) %>% 
  mutate(state=as.factor(state))

g <- tmp %>% 
  ggplot(aes(x=y1, y=y2)) +
  geom_point(data=tmp %>% filter(label=="truth"), colour="red") +
  geom_point(aes(shape=label), size=2) +
  xlab(latex2exp::TeX("$\\Y_{1t}")) +
  ylab(latex2exp::TeX("$\\Y_{2t}")) +
  scale_color_brewer("State", palette = "Dark2") +
  facet_wrap(~state, scales="free")
g

ggsave("../figures/state_cleanliness_mean.pdf", g)
Saving 5.64 x 3.48 in image

Output covariances for visualisation in Mathematica

Apply rubbish collection to real data

x <- read.csv("../data/hmm_test.csv", header = FALSE)

x <- x %>% 
  slice(1:nrow(x))
df_1 <- df %>% 
  filter(id == 10) %>% 
  mutate(dist = dist /max(dist)) %>% 
  mutate(time=seq_along(x)) %>% 
  slice(1:4000)

df_1 %>% 
  select(time, dist, angle) %>% 
  pivot_longer(-time) %>% 
  ggplot(aes(x=time, y=value)) +
  geom_line() +
  facet_wrap(~name, scales="free")

Fit normal model

data_stan <- list(
  N=nrow(df_1),
  dist=df_1$dist,
  angle=df_1$angle,
  K=2,
  cost=10,
  p_anomalous=0.01,
  conc_anomalous=100
)
model <- stan_model("hmm_wrapped_cauchy.stan")
DIAGNOSTIC(S) FROM PARSER:
Info:
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
    B_step[1] ~ normal(...)
Info:
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
    B_step[2] ~ normal(...)
fit <- optimise_repeat(5, model, data_stan)

df_1 %>% 
  select(time, dist, angle) %>% 
  mutate(state=as.factor(fit$par$state)) %>% 
  pivot_longer(-c(time, state)) %>% 
  ggplot(aes(x=time, y=value, colour=state)) +
  geom_line(aes(group=1)) +
  facet_wrap(~name, scales="free")

Using rubbish collector

data_stan <- list(
  N=nrow(df_1),
  dist=df_1$dist,
  angle=df_1$angle,
  K=2,
  cost=0.5,
  p_anomalous=0.1,
  conc_anomalous=10
)
model <- stan_model("hmm_bernoulli_wrapped_cauchy.stan")
DIAGNOSTIC(S) FROM PARSER:
Info:
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
    B_step[1] ~ normal(...)
Info:
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
    B_step[2] ~ normal(...)
fit <- optimise_repeat(2, model, data_stan)

df_1 %>% 
  select(time, dist, angle) %>% 
  mutate(state=as.factor(fit$par$state)) %>% 
  pivot_longer(-c(time, state)) %>% 
  ggplot(aes(x=time, y=value, colour=state)) +
  geom_line(aes(group=1)) +
  scale_color_brewer("State", palette = "Dark2") +
  facet_wrap(~name, scales="free")

g <- df_1 %>% 
  select(time, dist, angle) %>% 
  mutate(state=as.factor(fit$par$state)) %>% 
  pivot_longer(-c(time, state)) %>% 
  ggplot(aes(x=value, fill=state)) +
  stat_density(position="identity", alpha=0.5) +
  scale_fill_brewer("State", palette = "Dark2") +
  facet_wrap(~name, scales="free")
g
ggsave("../figures/angle_dist_separate_rubbish.pdf", g, width = 10, height = 4)

df_1 %>% 
  select(time, dist, angle) %>% 
  mutate(state=as.factor(fit$par$state)) %>% 
  ggplot(aes(x=dist, y=angle)) +
  geom_density_2d(aes(color=state))

LS0tCnRpdGxlOiAiSE1NIEJlcm5vdWxsaSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CnNvdXJjZSgiaGVscGVyLlIiKQpgYGAKCkdlbmVyYXRlIGRhdGEgd2l0aCBpbXBlcmZlY3Rpb25zCmBgYHtyfQptdXMgPC0gYygxLCA4KQpzaWdtYXMgPC0gYygxLCAxKQpkZnMgPC0gYygxLjIsIDEuMikKSyA8LSAyCm1fdHJhbnNpdGlvbl9wcm9icyA8LSBtYXRyaXgoYygwLjksIDAuMSwgMC4xLCAwLjkpLCBuY29sID0gMikKCm5fc3RlcHMgPC0gMTAwMApkZiA8LSBzaW11bGF0ZV9obW0obl9zdGVwcywgMSwgbV90cmFuc2l0aW9uX3Byb2JzLCAwLjIsIG11cywgc2lnbWFzLCBkZnMpCgpnIDwtIGRmICU+JSAKICBnZ3Bsb3QoYWVzKHg9dGltZSwgeT1vYnMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG91cj1zdGF0ZSkpICsKICB4bGFiKCJUaW1lIikgKwogIHlsYWIobGF0ZXgyZXhwOjpUZVgoIiRcXFlfdCIpKSArCiAgc2NhbGVfY29sb3JfYnJld2VyKCJTdGF0ZSIsIHBhbGV0dGUgPSAiRGFyazIiKQpnCgpnZ3NhdmUoIi4uL2ZpZ3VyZXMvc2ltdWxhdGVkX2NvcnJ1cHQucGRmIiwgZywgd2lkdGggPSA2LAogICAgICAgaGVpZ2h0ID0gMykKYGBgCgpGaXQgbm9ybWFsIEhNTSB3aXRoIDIgc3RhdGVzCmBgYHtyfQptb2RlbCA8LSBzdGFuX21vZGVsKCJobW0uc3RhbiIpCgpkYXRhX3N0YW4gPC0gbGlzdCgKICBOPW5yb3coZGYpLAogIGRpc3Q9ZGYkb2JzLAogIEs9MiwKICBjb3N0PTEwLAogIHBfYW5vbWFsb3VzPTAuMDEsCiAgY29uY19hbm9tYWxvdXM9MTAwCikKCmZpdCA8LSBvcHRpbWlzZV9yZXBlYXQoNSwgbW9kZWwsIGRhdGFfc3RhbikKCiMgc29tZSBpc3N1ZSB3aXRoIGFsZ28gZm9yIHQ9MQpkZiRzdGF0ZV9lc3QgPC0gZml0JHBhciRzdGF0ZQpkZiRzdGF0ZV9lc3RbMV0gPC0gMQoKZyA8LSBkZiAlPiUKICBtdXRhdGUoc3RhdGVfZXN0PWFzLmZhY3RvcihzdGF0ZV9lc3QpKSAlPiUgCiAgZ2dwbG90KGFlcyh4PXRpbWUsIHk9b2JzKSkgKwogIGdlb21fbGluZSgpICsKICBnZW9tX3BvaW50KGFlcyhjb2xvdXI9c3RhdGVfZXN0KSkgKwogIHhsYWIoIlRpbWUiKSArCiAgeWxhYihsYXRleDJleHA6OlRlWCgiJFxcWV90IikpICsKICBzY2FsZV9jb2xvcl9icmV3ZXIoIlN0YXRlIiwgcGFsZXR0ZSA9ICJEYXJrMiIpCgpnCmdnc2F2ZSgiLi4vZmlndXJlcy9zaW11bGF0ZWRfY29ycnVwdF9ub3JtYWxfc3RhdGUucGRmIiwgZywgd2lkdGggPSA2LAogICAgICAgaGVpZ2h0ID0gMykKYGBgCgpUcnkgbm9ybWFsIG1vZGVsIHdpdGggMyBzdGF0ZXMKYGBge3J9Cm1vZGVsIDwtIHN0YW5fbW9kZWwoImhtbS5zdGFuIikKCmRhdGFfc3RhbiA8LSBsaXN0KAogIE49bnJvdyhkZiksCiAgZGlzdD1kZiRvYnMsCiAgSz0zLAogIGNvc3Q9MTAsCiAgcF9hbm9tYWxvdXM9MC4wMSwKICBjb25jX2Fub21hbG91cz0xMDAKKQoKZml0IDwtIG9wdGltaXNlX3JlcGVhdCg1LCBtb2RlbCwgZGF0YV9zdGFuKQoKZGYkc3RhdGVfZXN0IDwtIGZpdCRwYXIkc3RhdGUKZGYkc3RhdGVfZXN0WzFdIDwtIDEKCmcgPC0gZGYgJT4lCiAgbXV0YXRlKHN0YXRlX2VzdD1hcy5mYWN0b3Ioc3RhdGVfZXN0KSkgJT4lIAogIGdncGxvdChhZXMoeD10aW1lLCB5PW9icykpICsKICBnZW9tX2xpbmUoKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPXN0YXRlX2VzdCkpICsKICB4bGFiKCJUaW1lIikgKwogIHlsYWIobGF0ZXgyZXhwOjpUZVgoIiRcXFlfdCIpKSArCiAgc2NhbGVfY29sb3JfYnJld2VyKCJTdGF0ZSIsIHBhbGV0dGUgPSAiRGFyazIiKQoKZwoKZ2dzYXZlKCIuLi9maWd1cmVzL3NpbXVsYXRlZF9jb3JydXB0X25vcm1hbF9zdGF0ZV8zLnBkZiIsIGcsIHdpZHRoID0gNiwKICAgICAgIGhlaWdodCA9IDMpCmBgYAoKCkZpdCBhIG1vZGVsIHdpdGggcnViYmlzaCBjb2xsZWN0aW9uIHN0YXRlLgpgYGB7cn0KZGF0YV9zdGFuIDwtIGxpc3QoCiAgTj1ucm93KGRmKSwKICBkaXN0PWRmJG9icywKICBLPTIsCiAgY29zdD0xMCwKICBwX2Fub21hbG91cz0wLjAxLAogIGNvbmNfYW5vbWFsb3VzPTEwMAopCm1vZGVsIDwtIHN0YW5fbW9kZWwoImhtbV9iZXJub3VsbGkuc3RhbiIpCmZpdCA8LSBvcHRpbWlzZV9yZXBlYXQoNSwgbW9kZWwsIGRhdGFfc3RhbikKCmRmJHN0YXRlX2VzdCA8LSBmaXQkcGFyJHN0YXRlCmRmJHN0YXRlX2VzdFsxXSA8LSAxCgpnIDwtIGRmICU+JQogIG11dGF0ZShzdGF0ZV9lc3Q9YXMuZmFjdG9yKHN0YXRlX2VzdCkpICU+JSAKICBnZ3Bsb3QoYWVzKHg9dGltZSwgeT1vYnMpKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG91cj1zdGF0ZV9lc3QpKSArCiAgeGxhYigiVGltZSIpICsKICB5bGFiKGxhdGV4MmV4cDo6VGVYKCIkXFxZX3QiKSkgKwogIHNjYWxlX2NvbG9yX2JyZXdlcigiU3RhdGUiLCBwYWxldHRlID0gIkRhcmsyIikKCmcKCmdnc2F2ZSgiLi4vZmlndXJlcy9zaW11bGF0ZWRfY29ycnVwdF9ydWJiaXNoLnBkZiIsIGcsIHdpZHRoID0gNiwKICAgICAgIGhlaWdodCA9IDMpCmBgYAoKYGBge3J9CmZpdCRwYXIkcGhpCmZpdCRwYXIkbXUKZml0JHBhciRzaWdtYQpgYGAKCgojIDJEIEhNTQpgYGB7cn0Kc3RhdGVzIDwtIHNpbXVsYXRlX2htbV9zdGF0ZXMobl9zdGVwcywgMSwgbV90cmFuc2l0aW9uX3Byb2JzKQpsaWJyYXJ5KG12dG5vcm0pCnJtdnJub3JtMkQgPC0gZnVuY3Rpb24obiwgbXV4LCBtdXksIHNpZ21heCwgc2lnbWF5LCByaG8pewogIHJldHVybihybXZub3JtKG4sIGMobXV4LCBtdXkpLAogICAgICAgICAgICAgICAgIG1hdHJpeChjKHNpZ21heF4yLCBzaWdtYXggKiBzaWdtYXkgKiByaG8sCiAgICAgICAgICAgICAgICAgICAgICAgICAgc2lnbWF4ICogc2lnbWF5ICogcmhvLCBzaWdtYXleMiksCiAgICAgICAgICAgICAgICAgICAgICAgIG5jb2wgPSAyKSkpCn0Kcm12dDJEIDwtIGZ1bmN0aW9uKG4sIG11eCwgbXV5LCBzaWdtYXgsIHNpZ21heSwgcmhvLCBkZil7CiAgcmV0dXJuKHJtdnQobiwKICAgICAgICAgICAgICBtYXRyaXgoYyhzaWdtYXheMiwgc2lnbWF4ICogc2lnbWF5ICogcmhvLAogICAgICAgICAgICAgICAgICAgICAgIHNpZ21heCAqIHNpZ21heSAqIHJobywgc2lnbWF5XjIpLAogICAgICAgICAgICAgICAgICAgICBuY29sID0gMiksCiAgICAgICAgICAgICAgZGYsCiAgICAgICAgICAgICAgYyhtdXgsIG11eSkpKQp9Cgp4IDwtIG1hdHJpeChucm93ID0gbGVuZ3RoKHN0YXRlcyksIG5jb2w9MikKcDFfY29ycnVwdGlvbiA8LSAwLjA1CnAyX2NvcnJ1cHRpb24gPC0gMC4xCmZvcihpIGluIHNlcV9hbG9uZyhzdGF0ZXMpKSB7CiAgdSA8LSBydW5pZigxKQogIGlmKHN0YXRlc1tpXSA9PSAxKSB7CiAgICBpZih1IDwgcDFfY29ycnVwdGlvbikgeyAKICAgICAgeFtpLCBdIDwtIHJtdnJub3JtMkQoMSwgLTUsIC01LCAyLCAyLCAwKQogICAgfSBlbHNlIHsKICAgICAgeFtpLCBdIDwtIHJtdnJub3JtMkQoMSwgMCwgMCwgMSwgMSwgMCkKICAgIH0KICB9IGVsc2UgewogICAgaWYodSA8IHAyX2NvcnJ1cHRpb24pIHsgCiAgICAgIHhbaSwgXSA8LSBybXZybm9ybTJEKDEsIDQsIDQsIDEsIDEsIDAuOCkKICAgIH0gZWxzZSB7CiAgICAgIHhbaSwgXSA8LSBybXZybm9ybTJEKDEsIDgsIDgsIDEsIDEsIDApCiAgICB9CiAgfQp9Cgp4IDwtIHggJT4lIAogIGFzLmRhdGEuZnJhbWUoKQoKZyA8LSB4ICU+JSAKICBtdXRhdGUoc3RhdGU9c3RhdGVzKSAlPiUgCiAgZ2dwbG90KGFlcyh4PVYxLCB5PVYyKSkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG91cj1hcy5mYWN0b3Ioc3RhdGUpKSkgKwogIHhsYWIobGF0ZXgyZXhwOjpUZVgoIiRcXFlfezF0fSIpKSArCiAgeWxhYihsYXRleDJleHA6OlRlWCgiJFxcWV97MnR9IikpICsKICBzY2FsZV9jb2xvcl9icmV3ZXIoIlN0YXRlIiwgcGFsZXR0ZSA9ICJEYXJrMiIpCmcKCmdnc2F2ZSgiLi4vZmlndXJlcy9zaW11bGF0ZWRfMmRfY29ycnVwdC5wZGYiLCBnLCB3aWR0aCA9IDYsCiAgICAgICBoZWlnaHQgPSAzKQpgYGAKCkZpdCAyRCByb2J1c3QgSE1NCmBgYHtyfQpkYXRhX3N0YW4gPC0gbGlzdCgKICBOPW5yb3coZGYpLAogIGRpc3Q9eCAlPiUgYXMubWF0cml4KCksCiAgSz0yLAogIGNvc3Q9NywKICBwX2Fub21hbG91cz0wLjEsCiAgY29uY19hbm9tYWxvdXM9MTAwCikKbW9kZWwgPC0gc3Rhbl9tb2RlbCgiaG1tX2Jlcm5vdWxsaV8yZC5zdGFuIikKZml0IDwtIG9wdGltaXNlX3JlcGVhdCgxMCwgbW9kZWwsIGRhdGFfc3RhbikKZml0X3IgPC0gZml0CgpnIDwtIHggJT4lIAogIG11dGF0ZShzdGF0ZT1maXRfciRwYXIkc3RhdGUpICU+JSAKICBnZ3Bsb3QoYWVzKHg9VjEsIHk9VjIpKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPWFzLmZhY3RvcihzdGF0ZSkpKSArCiAgeGxhYihsYXRleDJleHA6OlRlWCgiJFxcWV97MXR9IikpICsKICB5bGFiKGxhdGV4MmV4cDo6VGVYKCIkXFxZX3sydH0iKSkgKwogIHNjYWxlX2NvbG9yX2JyZXdlcigiU3RhdGUiLCBwYWxldHRlID0gIkRhcmsyIikKZwoKZ2dzYXZlKCIuLi9maWd1cmVzL3NpbXVsYXRlZF8yZF9jb3JydXB0X3J1YmJpc2gucGRmIiwgZywgd2lkdGggPSA2LAogICAgICAgaGVpZ2h0ID0gMykKYGBgCgpgYGB7cn0KZml0JHBhciRtdQpmaXQkcGFyJHNpZ21hCmBgYAoKVHJ5IHN0YW5kYXJkIEhNTSB3aXRoIDIgc3RhdGVzCmBgYHtyfQpkYXRhX3N0YW4gPC0gbGlzdCgKICBOPW5yb3coZGYpLAogIGRpc3Q9eCAlPiUgYXMubWF0cml4KCksCiAgSz0yLAogIGNvc3Q9NywKICBwX2Fub21hbG91cz0wLjEsCiAgY29uY19hbm9tYWxvdXM9MTAwCikKCm1vZGVsIDwtIHN0YW5fbW9kZWwoImhtbV8yZC5zdGFuIikKZml0IDwtIG9wdGltaXNlX3JlcGVhdCgxMCwgbW9kZWwsIGRhdGFfc3RhbikKZml0XzIgPC0gZml0CgojIG1hbnVhbCBzdGF0ZSBjbGVhbmluZwpzdGF0ZV8xIDwtIGZpdF8yJHBhciRzdGF0ZQpzdGF0ZV9jIDwtIGlmX2Vsc2Uoc3RhdGVfMT09MSwgMiwgMSkKCmcgPC0geCAlPiUgCiAgbXV0YXRlKHN0YXRlPXN0YXRlX2MpICU+JSAKICBnZ3Bsb3QoYWVzKHg9VjEsIHk9VjIpKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3VyPWFzLmZhY3RvcihzdGF0ZV9jKSkpICsKICB4bGFiKGxhdGV4MmV4cDo6VGVYKCIkXFxZX3sxdH0iKSkgKwogIHlsYWIobGF0ZXgyZXhwOjpUZVgoIiRcXFlfezJ0fSIpKSArCiAgc2NhbGVfY29sb3JfYnJld2VyKCJTdGF0ZSIsIHBhbGV0dGUgPSAiRGFyazIiKQpnCgpnZ3NhdmUoIi4uL2ZpZ3VyZXMvc2ltdWxhdGVkXzJkX2NvcnJ1cHRfbm9ybWFsXzIucGRmIiwgZywgd2lkdGggPSA2LAogICAgICAgaGVpZ2h0ID0gMykKYGBgCgpgYGB7cn0KZml0JHBhciRtdQpmaXQkcGFyJHNpZ21hCmBgYAoKClRyeSBzdGFuZGFyZCBITU0gd2l0aCAzIHN0YXRlcwpgYGB7cn0KZGF0YV9zdGFuIDwtIGxpc3QoCiAgTj1ucm93KGRmKSwKICBkaXN0PXggJT4lIGFzLm1hdHJpeCgpLAogIEs9MywKICBjb3N0PTcsCiAgcF9hbm9tYWxvdXM9MC4xLAogIGNvbmNfYW5vbWFsb3VzPTEwMAopCgptb2RlbCA8LSBzdGFuX21vZGVsKCJobW1fMmQuc3RhbiIpCmZpdCA8LSBvcHRpbWlzZV9yZXBlYXQoMTAsIG1vZGVsLCBkYXRhX3N0YW4pCmZpdF8zIDwtIGZpdAoKIyBtYW51YWwgc3RhdGUgY2xlYW5pbmcKc3RhdGVfMSA8LSBmaXRfMyRwYXIkc3RhdGUKc3RhdGVfYyA8LSBpZl9lbHNlKHN0YXRlXzE9PTEsIGlmX2Vsc2Uoc3RhdGVfMSksIDEpCgpnIDwtIHggJT4lIAogIG11dGF0ZShzdGF0ZT1zdGF0ZV9jKSAlPiUgCiAgZ2dwbG90KGFlcyh4PVYxLCB5PVYyKSkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG91cj1hcy5mYWN0b3Ioc3RhdGVfMSkpKSArCiAgeGxhYihsYXRleDJleHA6OlRlWCgiJFxcWV97MXR9IikpICsKICB5bGFiKGxhdGV4MmV4cDo6VGVYKCIkXFxZX3sydH0iKSkgKwogIHNjYWxlX2NvbG9yX2JyZXdlcigiU3RhdGUiLCBwYWxldHRlID0gIkRhcmsyIikKZwoKZ2dzYXZlKCIuLi9maWd1cmVzL3NpbXVsYXRlZF8yZF9jb3JydXB0X25vcm1hbF8zLnBkZiIsIGcsIHdpZHRoID0gNiwKICAgICAgIGhlaWdodCA9IDMpCmBgYAoKYGBge3J9CmZpdCRwYXIkbXUKZml0JHBhciRzaWdtYQpgYGAKClBsb3Qgc3RhdGUgY2xlYW5saW5lc3MKYGBge3J9Cm11X3IgPC0gZml0X3IkcGFyJG11Cm11X3IgPC0gbXVfcltjKDIsIDEpLCBdCm11XzIgPC0gZml0XzIkcGFyJG11Cm11XzIgPC0gbXVfMltjKDIsIDEpLCBdCm11XzMgPC0gZml0XzMkcGFyJG11W2MoMSwgMyksIF0KbXVfdHJ1ZSA8LSBtYXRyaXgoYygwLCAwLCA4LCA4KSwgbmNvbD0yLCBieXJvdyA9IFRSVUUpCgp0bXAgPC0gdHJpYmJsZSgKICB+eTEsIH55MiwgfnN0YXRlLCB+bGFiZWwsCiAgbXVfdHJ1ZVsxLCAxXSwgbXVfdHJ1ZVsxLCAyXSwgMSwgInRydXRoIiwKICBtdV90cnVlWzIsIDFdLCBtdV90cnVlWzIsIDJdLCAyLCAidHJ1dGgiLAogIG11X3JbMSwgMV0sIG11X3JbMSwgMl0sIDEsICJydWJiaXNoIiwKICBtdV9yWzIsIDFdLCBtdV9yWzIsIDJdLCAyLCAicnViYmlzaCIsCiAgbXVfMlsxLCAxXSwgbXVfMlsxLCAyXSwgMSwgIjItc3RhdGUgbm9ybWFsIiwKICBtdV8yWzIsIDFdLCBtdV8yWzIsIDJdLCAyLCAiMi1zdGF0ZSBub3JtYWwiLAogIG11XzNbMSwgMV0sIG11XzNbMSwgMl0sIDEsICIzLXN0YXRlIG5vcm1hbCIsCiAgbXVfM1syLCAxXSwgbXVfM1syLCAyXSwgMiwgIjMtc3RhdGUgbm9ybWFsIgopICU+JSAKICBtdXRhdGUoc3RhdGU9YXMuZmFjdG9yKHN0YXRlKSkKCmcgPC0gdG1wICU+JSAKICBnZ3Bsb3QoYWVzKHg9eTEsIHk9eTIpKSArCiAgZ2VvbV9wb2ludChkYXRhPXRtcCAlPiUgZmlsdGVyKGxhYmVsPT0idHJ1dGgiKSwgY29sb3VyPSJyZWQiKSArCiAgZ2VvbV9wb2ludChhZXMoc2hhcGU9bGFiZWwpLCBzaXplPTIpICsKICB4bGFiKGxhdGV4MmV4cDo6VGVYKCIkXFxZX3sxdH0iKSkgKwogIHlsYWIobGF0ZXgyZXhwOjpUZVgoIiRcXFlfezJ0fSIpKSArCiAgc2NhbGVfY29sb3JfYnJld2VyKCJTdGF0ZSIsIHBhbGV0dGUgPSAiRGFyazIiKSArCiAgZmFjZXRfd3JhcCh+c3RhdGUsIHNjYWxlcz0iZnJlZSIpCmcKCmdnc2F2ZSgiLi4vZmlndXJlcy9zdGF0ZV9jbGVhbmxpbmVzc19tZWFuLnBkZiIsIGcpCmBgYApPdXRwdXQgY292YXJpYW5jZXMgZm9yIHZpc3VhbGlzYXRpb24gaW4gTWF0aGVtYXRpY2EKYGBge3J9CnNpZ21hX3IgPC0gZml0X3IkcGFyJHNpZ21hCnNpZ21hXzIgPC0gZml0XzIkcGFyJHNpZ21hCnNpZ21hXzMgPC0gZml0XzMkcGFyJHNpZ21hCmBgYAoKCiMgQXBwbHkgcnViYmlzaCBjb2xsZWN0aW9uIHRvIHJlYWwgZGF0YQpgYGB7cn0KZGYgPC0gcmVhZC5jc3YoIi4uL2RhdGEvRWxlcGhhbnRzRGF0YV9hbm8uY3N2IixzZXA9IjsiKQpkZiA8LSBkZiAlPiUgCiAgcmVuYW1lX2FsbCh0b2xvd2VyKQpgYGAKCgpgYGB7cn0KZGZfMSA8LSBkZiAlPiUgCiAgZmlsdGVyKGlkID09IDEwKSAlPiUgCiAgbXV0YXRlKGRpc3QgPSBkaXN0IC9tYXgoZGlzdCkpICU+JSAKICBtdXRhdGUodGltZT1zZXFfYWxvbmcoeCkpICU+JSAKICBzbGljZSgxOjQwMDApCgpkZl8xICU+JSAKICBzZWxlY3QodGltZSwgZGlzdCwgYW5nbGUpICU+JSAKICBwaXZvdF9sb25nZXIoLXRpbWUpICU+JSAKICBnZ3Bsb3QoYWVzKHg9dGltZSwgeT12YWx1ZSkpICsKICBnZW9tX2xpbmUoKSArCiAgZmFjZXRfd3JhcCh+bmFtZSwgc2NhbGVzPSJmcmVlIikKYGBgCkZpdCBub3JtYWwgbW9kZWwKYGBge3J9CmRhdGFfc3RhbiA8LSBsaXN0KAogIE49bnJvdyhkZl8xKSwKICBkaXN0PWRmXzEkZGlzdCwKICBhbmdsZT1kZl8xJGFuZ2xlLAogIEs9MiwKICBjb3N0PTEwLAogIHBfYW5vbWFsb3VzPTAuMDEsCiAgY29uY19hbm9tYWxvdXM9MTAwCikKbW9kZWwgPC0gc3Rhbl9tb2RlbCgiaG1tX3dyYXBwZWRfY2F1Y2h5LnN0YW4iKQpmaXQgPC0gb3B0aW1pc2VfcmVwZWF0KDUsIG1vZGVsLCBkYXRhX3N0YW4pCgpkZl8xICU+JSAKICBzZWxlY3QodGltZSwgZGlzdCwgYW5nbGUpICU+JSAKICBtdXRhdGUoc3RhdGU9YXMuZmFjdG9yKGZpdCRwYXIkc3RhdGUpKSAlPiUgCiAgcGl2b3RfbG9uZ2VyKC1jKHRpbWUsIHN0YXRlKSkgJT4lIAogIGdncGxvdChhZXMoeD10aW1lLCB5PXZhbHVlLCBjb2xvdXI9c3RhdGUpKSArCiAgZ2VvbV9saW5lKGFlcyhncm91cD0xKSkgKwogIGZhY2V0X3dyYXAofm5hbWUsIHNjYWxlcz0iZnJlZSIpCmBgYAoKYGBge3J9CmRmXzEgJT4lIAogIHNlbGVjdCh0aW1lLCBkaXN0LCBhbmdsZSkgJT4lIAogIG11dGF0ZShzdGF0ZT1hcy5mYWN0b3IoZml0JHBhciRzdGF0ZSkpICU+JSAKICBwaXZvdF9sb25nZXIoLWModGltZSwgc3RhdGUpKSAlPiUgCiAgZ2dwbG90KGFlcyh4PXZhbHVlLCBmaWxsPXN0YXRlKSkgKwogIHN0YXRfZGVuc2l0eShwb3NpdGlvbj0iaWRlbnRpdHkiLCBhbHBoYT0wLjMpICsKICBmYWNldF93cmFwKH5uYW1lLCBzY2FsZXM9ImZyZWUiKQpgYGAKClVzaW5nIHJ1YmJpc2ggY29sbGVjdG9yCmBgYHtyfQpkYXRhX3N0YW4gPC0gbGlzdCgKICBOPW5yb3coZGZfMSksCiAgZGlzdD1kZl8xJGRpc3QsCiAgYW5nbGU9ZGZfMSRhbmdsZSwKICBLPTIsCiAgY29zdD0wLjUsCiAgcF9hbm9tYWxvdXM9MC4xLAogIGNvbmNfYW5vbWFsb3VzPTEwCikKbW9kZWwgPC0gc3Rhbl9tb2RlbCgiaG1tX2Jlcm5vdWxsaV93cmFwcGVkX2NhdWNoeS5zdGFuIikKZml0IDwtIG9wdGltaXNlX3JlcGVhdCgyLCBtb2RlbCwgZGF0YV9zdGFuKQoKZGZfMSAlPiUgCiAgc2VsZWN0KHRpbWUsIGRpc3QsIGFuZ2xlKSAlPiUgCiAgbXV0YXRlKHN0YXRlPWFzLmZhY3RvcihmaXQkcGFyJHN0YXRlKSkgJT4lIAogIHBpdm90X2xvbmdlcigtYyh0aW1lLCBzdGF0ZSkpICU+JSAKICBnZ3Bsb3QoYWVzKHg9dGltZSwgeT12YWx1ZSwgY29sb3VyPXN0YXRlKSkgKwogIGdlb21fbGluZShhZXMoZ3JvdXA9MSkpICsKICBzY2FsZV9jb2xvcl9icmV3ZXIoIlN0YXRlIiwgcGFsZXR0ZSA9ICJEYXJrMiIpICsKICBmYWNldF93cmFwKH5uYW1lLCBzY2FsZXM9ImZyZWUiKQpgYGAKCgpgYGB7cn0KZyA8LSBkZl8xICU+JSAKICBzZWxlY3QodGltZSwgZGlzdCwgYW5nbGUpICU+JSAKICBtdXRhdGUoc3RhdGU9YXMuZmFjdG9yKGZpdCRwYXIkc3RhdGUpKSAlPiUgCiAgcGl2b3RfbG9uZ2VyKC1jKHRpbWUsIHN0YXRlKSkgJT4lIAogIGdncGxvdChhZXMoeD12YWx1ZSwgZmlsbD1zdGF0ZSkpICsKICBzdGF0X2RlbnNpdHkocG9zaXRpb249ImlkZW50aXR5IiwgYWxwaGE9MC41KSArCiAgc2NhbGVfZmlsbF9icmV3ZXIoIlN0YXRlIiwgcGFsZXR0ZSA9ICJEYXJrMiIpICsKICBmYWNldF93cmFwKH5uYW1lLCBzY2FsZXM9ImZyZWUiKQpnCmdnc2F2ZSgiLi4vZmlndXJlcy9hbmdsZV9kaXN0X3NlcGFyYXRlX3J1YmJpc2gucGRmIiwgZywgd2lkdGggPSAxMCwgaGVpZ2h0ID0gNCkKYGBgCgoKYGBge3J9CmcgPC0gZGZfMSAlPiUgCiAgc2VsZWN0KHRpbWUsIGRpc3QsIGFuZ2xlKSAlPiUgCiAgbXV0YXRlKHN0YXRlPWFzLmZhY3RvcihmaXQkcGFyJHN0YXRlKSkgJT4lIAogIGdncGxvdChhZXMoeD1kaXN0LCB5PWFuZ2xlKSkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG91cj1zdGF0ZSksIHNpemU9MSkgKwogIHNjYWxlX2NvbG91cl9icmV3ZXIoIlN0YXRlIiwgcGFsZXR0ZSA9ICJEYXJrMiIpICsKICB4bGFiKCJEaXN0YW5jZSIpICsKICB5bGFiKCJBbmdsZSIpICsKICBmYWNldF93cmFwKH5zdGF0ZSkKZwpnZ3NhdmUoIi4uL2ZpZ3VyZXMvYW5nbGVfZGlzdF8yZF9ydWJiaXNoLnBkZiIsIGcsIHdpZHRoID0gMTAsIGhlaWdodCA9IDQpCmBgYAoKYGBge3J9CmRmXzEgJT4lIAogIHNlbGVjdCh0aW1lLCBkaXN0LCBhbmdsZSkgJT4lIAogIG11dGF0ZShzdGF0ZT1hcy5mYWN0b3IoZml0JHBhciRzdGF0ZSkpICU+JSAKICBnZ3Bsb3QoYWVzKHg9ZGlzdCwgeT1hbmdsZSkpICsKICBnZW9tX2RlbnNpdHlfMmQoYWVzKGNvbG9yPXN0YXRlKSkKYGBgCgo=